{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from pprint import pprint\n",
    "from skforecast.datasets import fetch_dataset\n",
    "from skforecast.drift_detection import PopulationDriftDetector, PopulationDriftDetector\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train: (8772, 11), Test: (8772, 11)\n"
     ]
    }
   ],
   "source": [
    "data = fetch_dataset('bike_sharing', verbose=False)\n",
    "data.head(3)\n",
    "\n",
    "data_train = data.iloc[: len(data)//2].copy()\n",
    "data_new  = data.iloc[len(data)//2 :].copy()\n",
    "print(f'Train: {data_train.shape}, Test: {data_new.shape}')\n",
    "data_train['weather'] = data_train['weather'].astype('category')\n",
    "data_new['weather'] = pd.Categorical(data_new['weather'], categories=data_train['weather'].cat.categories)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\jaesc2\\AppData\\Local\\Temp\\ipykernel_39756\\2547271043.py:57: UserWarning: Dataset has 0 variance; skipping density estimate. Pass `warn_singular=False` to disable this warning.\n",
      "  sns.kdeplot(distances, ax=ax_hist, fill=True)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GIF saved as population_drift_detection.gif\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "Could not find a backend to open `population_drift_detection.mp4`` with iomode `wI`.\nBased on the extension, the following plugins might add capable backends:\n  FFMPEG:  pip install imageio[ffmpeg]\n  pyav:  pip install imageio[pyav]",
     "output_type": "error",
     "traceback": [
      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
      "\u001b[31mValueError\u001b[39m                                Traceback (most recent call last)",
      "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[5]\u001b[39m\u001b[32m, line 111\u001b[39m\n\u001b[32m    109\u001b[39m \u001b[38;5;66;03m# Save as MP4 video\u001b[39;00m\n\u001b[32m    110\u001b[39m mp4_filename = \u001b[33m'\u001b[39m\u001b[33mpopulation_drift_detection.mp4\u001b[39m\u001b[33m'\u001b[39m\n\u001b[32m--> \u001b[39m\u001b[32m111\u001b[39m \u001b[43mimageio\u001b[49m\u001b[43m.\u001b[49m\u001b[43mmimsave\u001b[49m\u001b[43m(\u001b[49m\n\u001b[32m    112\u001b[39m \u001b[43m    \u001b[49m\u001b[43mmp4_filename\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m    113\u001b[39m \u001b[43m    \u001b[49m\u001b[43mframes\u001b[49m\u001b[43m,\u001b[49m\n\u001b[32m    114\u001b[39m \u001b[43m    \u001b[49m\u001b[43mfps\u001b[49m\u001b[43m=\u001b[49m\u001b[43mfps\u001b[49m\u001b[43m,\u001b[49m\u001b[43m  \u001b[49m\u001b[38;5;66;43;03m# frame rate\u001b[39;49;00m\n\u001b[32m    115\u001b[39m \u001b[43m    \u001b[49m\u001b[43mcodec\u001b[49m\u001b[43m=\u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43mlibx264\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m  \u001b[49m\u001b[38;5;66;43;03m# high-quality H.264 codec\u001b[39;49;00m\n\u001b[32m    116\u001b[39m \u001b[43m    \u001b[49m\u001b[43mquality\u001b[49m\u001b[43m=\u001b[49m\u001b[32;43m10\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m       \u001b[49m\u001b[38;5;66;43;03m# 10 = best quality, 0 = worst\u001b[39;49;00m\n\u001b[32m    117\u001b[39m \u001b[43m    \u001b[49m\u001b[43mffmpeg_params\u001b[49m\u001b[43m=\u001b[49m\u001b[43m[\u001b[49m\n\u001b[32m    118\u001b[39m \u001b[43m        \u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43m-crf\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43m17\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m         \u001b[49m\u001b[38;5;66;43;03m# lower CRF = higher quality (range: 0–51)\u001b[39;49;00m\n\u001b[32m    119\u001b[39m \u001b[43m        \u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43m-pix_fmt\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43myuv420p\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[38;5;66;43;03m# ensures wide compatibility\u001b[39;49;00m\n\u001b[32m    120\u001b[39m \u001b[43m        \u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43m-preset\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[33;43m'\u001b[39;49m\u001b[33;43mslow\u001b[39;49m\u001b[33;43m'\u001b[39;49m\u001b[43m     \u001b[49m\u001b[38;5;66;43;03m# better compression efficiency (options: ultrafast → placebo)\u001b[39;49;00m\n\u001b[32m    121\u001b[39m \u001b[43m    \u001b[49m\u001b[43m]\u001b[49m\n\u001b[32m    122\u001b[39m \u001b[43m)\u001b[49m\n\u001b[32m    123\u001b[39m \u001b[38;5;28mprint\u001b[39m(\u001b[33mf\u001b[39m\u001b[33m'\u001b[39m\u001b[33mHigh-quality MP4 saved as \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mmp4_filename\u001b[38;5;132;01m}\u001b[39;00m\u001b[33m'\u001b[39m)\n",
      "\u001b[36mFile \u001b[39m\u001b[32mc:\\Users\\jaesc2\\Miniconda3\\envs\\skforecast_py12\\Lib\\site-packages\\imageio\\v2.py:494\u001b[39m, in \u001b[36mmimwrite\u001b[39m\u001b[34m(uri, ims, format, **kwargs)\u001b[39m\n\u001b[32m    492\u001b[39m imopen_args = decypher_format_arg(\u001b[38;5;28mformat\u001b[39m)\n\u001b[32m    493\u001b[39m imopen_args[\u001b[33m\"\u001b[39m\u001b[33mlegacy_mode\u001b[39m\u001b[33m\"\u001b[39m] = \u001b[38;5;28;01mTrue\u001b[39;00m\n\u001b[32m--> \u001b[39m\u001b[32m494\u001b[39m \u001b[38;5;28;01mwith\u001b[39;00m \u001b[43mimopen\u001b[49m\u001b[43m(\u001b[49m\u001b[43muri\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[33;43m\"\u001b[39;49m\u001b[33;43mwI\u001b[39;49m\u001b[33;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m*\u001b[49m\u001b[43m*\u001b[49m\u001b[43mimopen_args\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mas\u001b[39;00m file:\n\u001b[32m    495\u001b[39m     \u001b[38;5;28;01mreturn\u001b[39;00m file.write(ims, is_batch=\u001b[38;5;28;01mTrue\u001b[39;00m, **kwargs)\n",
      "\u001b[36mFile \u001b[39m\u001b[32mc:\\Users\\jaesc2\\Miniconda3\\envs\\skforecast_py12\\Lib\\site-packages\\imageio\\core\\imopen.py:281\u001b[39m, in \u001b[36mimopen\u001b[39m\u001b[34m(uri, io_mode, plugin, extension, format_hint, legacy_mode, **kwargs)\u001b[39m\n\u001b[32m    275\u001b[39m         err_msg += (\n\u001b[32m    276\u001b[39m             \u001b[33m\"\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[33mBased on the extension, the following plugins might add capable backends:\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[33m\"\u001b[39m\n\u001b[32m    277\u001b[39m             \u001b[33mf\u001b[39m\u001b[33m\"\u001b[39m\u001b[38;5;132;01m{\u001b[39;00minstall_candidates\u001b[38;5;132;01m}\u001b[39;00m\u001b[33m\"\u001b[39m\n\u001b[32m    278\u001b[39m         )\n\u001b[32m    280\u001b[39m request.finish()\n\u001b[32m--> \u001b[39m\u001b[32m281\u001b[39m \u001b[38;5;28;01mraise\u001b[39;00m err_type(err_msg)\n",
      "\u001b[31mValueError\u001b[39m: Could not find a backend to open `population_drift_detection.mp4`` with iomode `wI`.\nBased on the extension, the following plugins might add capable backends:\n  FFMPEG:  pip install imageio[ffmpeg]\n  pyav:  pip install imageio[pyav]"
     ]
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "from skforecast.plot import set_dark_theme\n",
    "set_dark_theme()\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "import os\n",
    "from PIL import Image\n",
    "import numpy as np\n",
    "\n",
    "try:\n",
    "    import imageio\n",
    "except ImportError:\n",
    "    print(\"imageio not installed. Install with: pip install imageio\")\n",
    "    raise\n",
    "\n",
    "feature = 'temp'\n",
    "ref_data = data_train[feature]\n",
    "distances = []\n",
    "fps = 1  # frames per second for video\n",
    "output_width = 900  # width for resized frames\n",
    "frames = []\n",
    "\n",
    "# Generate chunk starts every 2 months starting from January 2011\n",
    "chunk_starts = pd.date_range(start='2011-01-01', end=data_train.index.max(), freq='1MS')\n",
    "\n",
    "for i, chunk_start in enumerate(chunk_starts[:-1]):  # Exclude the last one to avoid going beyond data\n",
    "    chunk_end = min(chunk_start + pd.DateOffset(months=1) - pd.Timedelta(hours=1), data_train.index.max())\n",
    "    data_chunk = data_train.loc[chunk_start:chunk_end, feature]\n",
    "    \n",
    "    # Plot time series (spanning first row)\n",
    "    fig = plt.figure(figsize=(12, 8))\n",
    "    ax_ts = plt.subplot2grid((2, 2), (0, 0), colspan=2)\n",
    "    data_train.loc[:, 'temp'].plot(ax=ax_ts, label='Reference data')\n",
    "    data_chunk.plot(ax=ax_ts, label=f'Chunk {i}', color='red')\n",
    "    ax_ts.set_xlabel('')\n",
    "    ax_ts.set_ylabel('')\n",
    "    ax_ts.set_title('Distance-Based framework for temporal drift detection',\n",
    "                     fontsize=18, pad=25, fontweight='semibold')\n",
    "    ax_ts.legend(loc='upper left')\n",
    "    \n",
    "    # Compare distributions using kdeplot (second row, first column)\n",
    "    ax_kde = plt.subplot2grid((2, 2), (1, 0))\n",
    "    chunk_data = data_new.loc[chunk_start:chunk_end, feature]\n",
    "    sns.kdeplot(ref_data, label='Reference data', color='#30a2da', ax=ax_kde)\n",
    "    sns.kdeplot(data_chunk, label=f'Chunk {i}', color='red', ax=ax_kde)\n",
    "    ax_kde.set_title(f'Distribution comparison for chunk {i}')\n",
    "    ax_kde.set_xlabel('')\n",
    "    #ax_kde.legend(loc='upper left')\n",
    "\n",
    "    # Histogram of calculated distances (second row, second column)\n",
    "    from scipy.stats import ks_2samp\n",
    "    ks_statistic, p_value = ks_2samp(ref_data, data_chunk)\n",
    "    distances.append(ks_statistic)\n",
    "    \n",
    "    # Empirical distribution of distances\n",
    "    ax_hist = plt.subplot2grid((2, 2), (1, 1))\n",
    "    sns.kdeplot(distances, ax=ax_hist, fill=True)\n",
    "    # add rug plot\n",
    "    sns.rugplot(distances, ax=ax_hist, color='white')\n",
    "    ax_hist.set_title('Empirical distribution of distances')\n",
    "    ax_hist.set_xlabel('')\n",
    "    ax_hist.set_ylabel('')\n",
    "\n",
    "    # If it is the last chunk, plot final histogram with highlighted quantiles\n",
    "    if i == len(chunk_starts) - 2:\n",
    "        quantile_95 = np.quantile(distances, 0.95)\n",
    "        ax_hist.axvline(quantile_95, color='white', linestyle='--', label='95th Percentile')\n",
    "\n",
    "        # Add annotation with arrow\n",
    "        ax_hist.annotate(\n",
    "            '95th Percentile',\n",
    "            xy=(quantile_95, ax_hist.get_ylim()[1] * 0.8),       # Point to the line\n",
    "            xytext=(quantile_95 + (ax_hist.get_xlim()[1] - ax_hist.get_xlim()[0]) * 0.05, \n",
    "                    ax_hist.get_ylim()[1] * 0.9),                 # Place text slightly to the right/top\n",
    "            arrowprops=dict(facecolor='white', shrink=0.05, width=1, headwidth=6),\n",
    "            color='white',\n",
    "            fontsize=12,\n",
    "            ha='left',\n",
    "            va='center'\n",
    "        )\n",
    "    \n",
    "    fig.tight_layout()\n",
    "    #plt.show();\n",
    "\n",
    "    # Save frame\n",
    "    filename = f\"temp_{i}.png\"\n",
    "    fig.savefig(filename, dpi=120, bbox_inches='tight')\n",
    "    plt.close(fig)\n",
    "\n",
    "    img = Image.open(filename)\n",
    "    img = img.resize((output_width, int(output_width * img.height / img.width)))\n",
    "    frames.append(np.array(img))\n",
    "    os.remove(filename)\n",
    "\n",
    "# Save as GIF (slower and loop infinitely)\n",
    "gif_filename = 'population_drift_detection.gif'\n",
    "# Duplicate the last frame to extend its display time\n",
    "frames.extend([frames[-1]] * 5)\n",
    "\n",
    "imageio.mimsave(\n",
    "    gif_filename,\n",
    "    frames,\n",
    "    format='GIF',\n",
    "    duration=900,  # 900 milliseconds per frame\n",
    "    loop=0\n",
    ")\n",
    "print(f'GIF saved as {gif_filename}')\n",
    "\n",
    "# Save as MP4 video\n",
    "mp4_filename = 'population_drift_detection.mp4'\n",
    "imageio.mimsave(\n",
    "    mp4_filename,\n",
    "    frames,\n",
    "    fps=fps,  # frame rate\n",
    "    codec='libx264',  # high-quality H.264 codec\n",
    "    quality=10,       # 10 = best quality, 0 = worst\n",
    "    ffmpeg_params=[\n",
    "        '-crf', '17',         # lower CRF = higher quality (range: 0–51)\n",
    "        '-pix_fmt', 'yuv420p',# ensures wide compatibility\n",
    "        '-preset', 'slow'     # better compression efficiency (options: ultrafast → placebo)\n",
    "    ]\n",
    ")\n",
    "print(f'High-quality MP4 saved as {mp4_filename}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "skforecast_py12",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.11"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
